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PREFACE 


The  work  reported  herein  was  supported  by  the  US  Army  Ballistic 
Research  Laboratory  (BRL),  United  States  Army  Research  and  Development  Command 
(ARRADCOM)  and  performed  by  Science  Applications,  Inc.  (SAI)  under  Contract 
DAAK11-80-C-0079  dated  30  Jul  80.  The  contract  was  supported  in  part  by  funds 
provided  by  the  United  States  Air  Force  Weapons  Laboratory  (AFWL).  The  AFWL 
is  interested  in  ultimately  making  this  new  approach  available  in  the  Veotor 
HULL  code,  operational  on  the  AFWL  Cray-1.  The  BRL  version  of  HULL  is  opera¬ 
tional  on  a  Control  Data  Corporation  (CDC)  7600  computer  at  BRL. 

The  work  performed  under  this  contract  consisted  of:  development 
of  computer  algorithms  to  treat  partially  reflective  cells  in  two  dimensions 
(both  cylindrical  and  Cartesian),  implementation  of  these  algorithms  in  the 
BRL  Terminal  Ballistics  Division  (TBD)  version  of  the  AFWL  HULL  code,  and 
support  to  BRL  in  debugging  of  the  implementation  and  test  of  the  algorithms. 

Mr.  Richard  E.  Lottero  was  the  BRL  contracting  officer's 
technical  representative.  SAI  personnel  worked  closely  with  him  and 
Mr.  John  D.  Wortman,  also  of  BRL.  Mr.  Burton  S.  Chambers  III  of  SAI  was  the 
principal  investigator  for  this  effort.  Dr.  John  A.  Hasdai  assisted  by 
providing  guidance  on  how  to  implement  these  algorithms  for  future  compati¬ 
bility  with  the  Vector  Hull  code  being  developed  for  the  AFWL  by  SAI  under 
another  contract. 

The  first  author  would  like  to  extend  his  sincere  appreciation  to 
Mr.  Wortman  for  his  helpful  suggestions  and  Improvements  to  some  of  the 
algorithms.  His  help  was  essential  to  the  successful  completion  of  this 
effort,  since  he  performed  most  of  the  check-out  calculations.  The  authors 
extend  their  slncerest  thanks  to  Mr.  Lottero  and  Dr.  Clarence  W.  Kitchens,  Jr. 
for  continuing  guidance  and  interest,  without  which  this  effort  would  not  have 
been  possible. 
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Section  1 
INTRODUCTION 

In  order  to  allow  realistic  hardening  design  criteria  to  be  estab¬ 
lished  for  the  ultimate  protection  and  survivability  of  U.S.  Army  military 
equipment  (e.g.,  communications  shelters  and  antennae,  from  nuclear  warfare 
environments,  primarily  blast  effects),  it  is  necessary  to  have  reasonably 
accurate  estimates  of  the  time-history  of  the  airblast-induced  loads. 
Although  much  experimental  work  has  been  done  to  produce  high-fidelity  simula¬ 
tions  of  these  loads,  there  exists  a  need  for  the  ability  to  predict  these 
loads  theoretically.  Current  simple  predictive  models  are  at  times  inade¬ 
quate,  and  more  detailed  flow  field  interaction  calculations  are  advisable. 
Unfortunately,  the  U.S.  Army  has  been  hampered  in  its  ability  to  predict 
detailed  blast  loading  pressure  distributions  on  such  structures  because  their 
computational  tools  do  not  adequately  treat  the  complex  geometry  of  these 
structures . 


On  the  scale  of  the  nuclear  event,  it  is  a  reasonable  approxima¬ 
tion  to  treat  the  airblast  flow  as  inviscid,  and  treat  objects  in  the  flow 
field's  vicinity  as  rigid  reflective  bodies  during  the  diffraction  phase  of 
the  shock  loading.  This  is  the  approach  historically  taken  by  BRL,  AFWL,  and 
SAI  for  predicting  early-time  airblast  loading  on  structures.  The  results  for 
many  cases  have  been  encouraging  and  in  some  cases  extremely  helpful;  however, 
room  for  improvement  exists. 

The  experience  at  BRL  with  the  HULL  code  (references  1  through  4), 
and  independently  at  SAI  (references  5  and  6),  has  demonstrated  that  HULL,  an 
efficient  multi-material,  multi-dimensional  inviscid  hydrodynamics  code  with 
an  option  for  treating  structures  as  perfectly  reflective  cells,  can  success¬ 
fully  predict  blast  loading  on  generic  non-responding  shapes  if  the  target 
surfaces  conform  with  flow  field  cell  boundaries.  It  was  also  recognized  at 
both  BRL  and  SAI  that  there  existed  a  need  to  modify  HULL  to  improve  its  capa¬ 
bility  for  irregular  rigid  structures  where  the  above  condition  was  not  met. 

This  effort,  therefore,  was  performed  to  improve  the  BRL  HULL 

code's  ability  for  predicting  loads  on  surfaces  not  parallel  with  the 
coordinate  axes.  The  approach  is  to  incorporate  into  HULL  a  capability  for 
treating  cells  which  are  partially  fluid  and  partially  reflective.  The  effort 
was  to  consist  of  the  following  tasks:  (1)  Development  of  a  numerical 

algorithm  for  treating  two-dimensional  partial  hydrodynamic/partial  rigid 

cells,  (2)  Implementation  of  this  algorithm  in  the  2-D  version  (Cartesian  and 
cylindrical)  of  BRL  HULL,  and  (3)  Implementation  of  a  similar  algorithm  in  the 
3-D  version  of  BRL  HULL.  Because  of  the  complexity  of  the  effort  and  the 
uncertainty  associated  with  specifying  operational  modifications  to  2-D  BRL 
HULL,  Task  3  was  necessarily  predicated  on  finishing  Tasks  1  and  2. 

The  first  task  consisted  of  the  design  and  development  of  numer¬ 
ical  algorithms  compatible  with  the  BRL  HULL  computer  program  for  treating 
hybrid  (hereafter  referred  to  as  shore)  computational  cells  that  are  partly 
hydrodynamic  and  partly  rigid  material.  This  extended  the  "island"  concept 
presently  in  HULL.  The  techniques  developed  treat  the  reflective  boundary 

conditions  in  the  shore  cells,  and  update  the  computation  of  the  flow  vari¬ 
ables  in,  and  in  the  vicinity  of,  tnese  shore  cells.  The  algorithms  were 
formulated  in  both  Cartesian  and  cylindrical  geometry  for  implementation  in 


the  2-D  version  of  BRL  HULL.  Consideration  was  given  to  any  restrictions 
imposed  on  the  fluxing  algorithms  and  the  time-step  algorithms  by  these  shore 
cells. 

The  second  task  consisted  of  the  design  of  code  architectural  mod¬ 
ifications  necessary  to  implement  the  treatment  of  shore  cells  in  the  2-D 
version  of  BRL  HULL.  Various  techniques  for  distinguishing  shore  cells  were 
considered,  and  the  method  selected  was  considered  the  best  one  primarily 
because  of  its  ease  of  implementation,  and  efficient  calculation.  All  modifi¬ 
cations  were  implemented  in  the  form  of  a  HULL  option,  named  SHORE,  to  avoid 
degrading  the  operational  status  of  HULL  and  to  avoid  introducing  extraneous 
coding.  This  latter  design  goal  of  HULL  directly  contributes  to  its  speed  of 
operation.  Furthermore,  all  architectural  modifications  are  considered  to  be 
consistent  with  the  Vector  HULL  architecture,  and  hence,  are  consistent  with 
the  vector  architecture  of  a  Cray-1  vector  processor. 

Under  Task  2,  SAI  also  assisted  BRL  personnel  in  setting  up  and 
running  2-D  test  calculations.  These  calculations  were  designed  to  allow  a 
determination  to  be  made  of:  (1)  the  correctness  of  the  modifications  and  (2) 
the  agreement  with  some  available  experimental  data. 

The  third  task  was  to  have  consisted  of  implementing  a  modified 
version  of  the  2-D  Cartesian  SHORE  algorithm  into  3-D  HULL.  However,  Tasks  1 
and  2  required  more  effort  than  originally  planned,  with  several  variations  on 
the  algorithm  tried  before  a  satisfactory  one  was  found.  The  work  performed 
under  Task  3  did  indicate  that  the  implementation  of  a  modified  algorithm  for 
3-D  SHORE  cells  is  feasible,  but  will  require  more  effort  than  originally 
estimated. 
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Section  2 

CURRENT  HULL  CODE  FORMULATION 


This  section  describes  the  current  approach  in  the  BRL  HULL  code  to 
solving  the  partial  differential  equations  describing  inviscid  fluid  flow.  It 
also  describes  the  approach  taken  to  represent  rigid  structures  before  this 
effort  was  performed.  This  approach  was  to  introduce  cells  into  HULL  that  are 
non-hydrodynamic  and  perfectly  reflective.  These  reflective  cells  are  called 
islands,  and  an  ensemble  of  them  simulates  a  solid  structure. 

2.1  The  Differential  Equations 

The  HULL  code  is  designed  to  efficiently  solve  the  hyperbolic 
partial  differential  equations  describing  inviscid,  nonconducting  fluid  flow 
in  the  form: 


$ 


u  =  0 


(1) 


P  f  .  =  p| 


p  +  V  •  (pu)  =  p  u  *  g 
dt 


1  -Y  -*■ 

E  a  I  +  2  u  *  u 


(2) 

(3) 

(4) 


where 


P 

P 

■> 

u 

I 

g 

t 


3 

material  density  (g/cra  ) 

2 

pressure  (dynes/cm  ) 

(u,v,w)  the  fluid  velocity  (cm/s) 
specific  internal  energy  (ergs/g) 

p 

acceleration  of  gravity  (cm/s  ) 
time  (s). 


These  differential  equations  can  be  solved  using  various  numerical 
difference  schemes.  In  the  current  implementations  of  HULL,  the  difference 
equations  are  formulated  in  either  Cartesian  or  cylindrical  coordinates  for 
the  2-D  versions  and  in  Cartesian  coordinates  for  the  3-D  version.  The 
difference  equations  for  3-D  are  described  below.  They  represent  the  solution 
for  the  above  equations  when  the  independent  spatial  coordinates  are  so-called 
Lagrangian  coordinates,  that  is,  they  move  with  the  fluid.  The  difference 
method  is  an  explicit  conservative  method,  which  specifically  could  be  consid¬ 
ered  a  modification  of  a  two-step  Lax-Wendroff  (reference  7)  scheme.  Since 
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the  original  partial  differential  equations  are  classified  as  hyperbolic,  a 
Courant-Friedrichs-Lewy  condition  is  imposed  on  the  time  step  to  assure  numer¬ 
ical  stability. 

Although  the  finite  difference  analogs  to  the  above  equations  could 
be  written  down  directly,  another  equation  is  used  in  the  differencing.  This 
equation: 


d£ 

dt 


*  p(Vf’ 


( V»u) 


0, 


(5) 


where  Yeff  is  the  effective  ratio  of  specific  heats,  is  derived  in  Appendix  A. 


Although  HULL  has  been  used  quite  successfully  for  a  variety  of 
problems,  it  should  be  understood  that  while  equation  (5)  is  shown  to  be 
exact,  it  is  used  in  HULL  in  an  approximate  way  in  that  it  is  assumed  that  the 
effective  gamma  of  the  gas,  calculated  as 


eff  = 


(6) 


is  constant  over  the  interval  of  the  time  step.  This  approximation,  which 
amounts  to  ignoring  the  time  derivative  of  gamma,  is  not  strictly  valid  for 
all  regimes  and  warrants  additional  investigation.  Nevertheless,  this  approx¬ 
imation  has  been  historically  used  for  HULL  air  blast  calculations. 

The  remainder  of  this  section  consists  of  a  presentation  of  the 
difference  equations  typically  used  for  3-D  air  blast  problems.  The  formula¬ 
tion  given  is  as  found  in  the  BRL  HULL  code  and  is  included  as  a  prelude  to 
the  discussions  on  the  shore  cell  concept.  Approximations  made  in  solving  the 
original  differential  equations  are  also  identified.  Only  the  first  phase  of 
the  usual  HULL  technique  is  included  here.  The  first  phase  solves  the  differ¬ 
ential  equations  couched  in  a  Lagrangian  frame  of  reference,  that  is,  the 
frame  moves  with  the  fluid.  (Although  HULL  is  an  Eulerian  code,  the  approach 
is  to  calculate  in  the  Lagrangian  reference  frame,  and  then  to  move  mass, 
momentum,  and  energy  rather  than  material  boundaries.)  The  second  phase  uses  a 
donor  cell  technique. 

2.2  The  3-D  HULL  Differencing  Scheme 

The  differencing  in  the  Lagrangian  phase  is  a  two  step  technique. 
The  equations  are  presented  in  reverse  order  of  calculation. 

Equation  (3)  repeated  here  for  convenience,  is 


„dE  n 

pdt 


( pu )  =  p  U*  g 


In  a  3-D  Cartesian  frame  of  reference  the  divergence  is 


where  u,  v,  and  w  are  the  velocity  components  in  the  x,  y,  and  z  directions, 
respectively. 


‘  (PW)  i,j,k+J|“"iwj 


where  is  a  gravitational  potential  term  and  E  is  the  total  energy  (i.e., 
kinetic  +  internal).  The  mass  in  cell  (i,j,k)  at  time  (n)  is  denoted  by 
m?,J.k*  The  superscript  n  means  the  value  is  at  the  current  time,  and  n+1 
denotes  the  value  at  the  updated  time  (or  a  one  time  step  advance  from  n  to 
n+1).  Similarily  the  superscript  n+J  means  the  value  is  at  the  center  of  the 
time  step.  The  subscripts  indicate  positions  in  space;  the  integer  values 
denote  cell  centers,  differences  of  one  denote  adjacent  cells,  and  half 
integer  values  means  the  quantity  is  at  one  of  the  boundaries  of  the  cell. 
This  same  notation  is  used  even  when  the  cells  are  not  of  equal  size. 


The  cell  denoted  by  i,  j,  and  k  lies  between  x^_1  and 
y^_1  and  y^,  and  between  zk-1  and  zk-  Further,  fai  =  x^-x.^,  Ay ^  = 


and  fa, 


zk"zk-1* 


between 


Two  confusing  features  of  the  difference  equations  in  HULL  are 
illustrated  by  Equation  (7).  Consider 


]_  9(pu) 
P  3x 


n+  £ 

i-L  J.k 


-  (pu) 


n+  i 

i-*-L  jik 
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n 

i.  j.k 


Ay  .fa 
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The  negative  sign  is  supplied  by  a  reversal  of  the  natural  order  of  differ¬ 
encing  and 


Ay .  fa 

Jj _ * 


n 

m. 


Ax.p.  .  . 
i ,  j ,  k  1  i.J.k 


The  latter  is  true  because  m.  .  .  ,  the  mass  of  fluid  in  cell  (i,j,k)  at  time 

1  >  J  9  K  - 

(n),  is  equal  to  the  product  of  density,  P.  ,  .  and  volume,  V.  =  Ax  Ay  Az  . 

1 t J t K  1 JK  1  J  K 

The  velocities  are  updated  to  time  (n+1)  by  differencing  equation 
(2),  which  is 


_  du  n  * 
P  dt  +  Vp  =  pg 


This  is  differenced  for  each  scalar  component  of  velocity  as: 
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In  order  to  evaluate  these  expressions  the  values  for  the  half  time 
step  (e.g.,  pn+£  and  (pu)n+J)  are  needed  on  the  boundaries  of  the  unit  cell. 
The  pressure  at  time  (n+J)  is  obtained  from  the  differential  equation  derived 
in  Appendix  A  and  repeated  here: 

£  +  P  W.crKV  •  I)  =  0 


3?  *  'P  <W 


f9u  3v  3w"l 
|3x  +  3y  +  3zJ 


in  3-D  Cartesian,  coordinates. 


In  order  to  solve  this  equation  without  computing  hydrodynamic 
quantities  on  the  cell  corners,  an  approximation  is  made  in  HULL.  The  partial 
derivatives  at  a  given  side  (i.e.,  a  boundary)  of  a  cell  are  assumed  to  be 
small  in  the  plane  of  the  side  and  are  ignored.  While  this  approximation  has 
been  used  since  HULL  was  first  implemented,  we  believe  additional  work  is 
needed  in  understanding  the  effect  of  this  approximation.  HULL  differences 
the  above  equation  at  cell  boundary  i+|,  j,  k  as: 

n+i  n  _ At  /  n  \f  n  n  1)  (11) 

Pi+J,j,k  -  pW,j,k  -  2Axw  (  Ye  ff1+i|J|k)Ui+1,J,k  " 


where 


Ax.  , 
i+5 


AXi  +  Axi+1 


The  values  of  u  ,  ...  on  the  boundaries  are  obtained  in  a  manner 
similar  to  what  was  done  for  the  velocity  at  the  whole  step;  however,  here  the 
differencing  (assuming  equally-sized  cells)  is  occurring  at  the  boundaries: 


W,  J,k 


n 

UW,  j,k 


i-^1  -  pi,j,k^ 

Ax.  ,  P^+i  .  . 
i+5  i+5,  j,k 


-/here 


n+J  _  n 
Pi+J,j,k  “  Pi+i,i,k 


1.0  -  M 
2 


and 


n  n 

n  .  mi,j,k  *  mi+1 ,  j,k  (14) 

pi+$,j,k  “  (Axt  +  Axi+1 )  Ayj  Azk* 

The  product  (pu)n+^  in  Equation  (7)  is  obtained  directly  by 

(pu)n+*  =  pn+*un+*.  (15) 

and  (pv)n+^  and  (pw)n+^  are  similarly  computed. 


In  order  to  conserve  storage,  reduce  data  handling,  and  increase 
efficiency,  an  orderly  sweep  is  made  through  the  mesh,  advancing  the  cell 
hydrodynamic  quantities  at  the  appropriate  time  when  the  necessary  data  are 
available,  but  being  careful  not  to  destroy  a  quantity  if  it  is  needed  later. 
To  accomplish  this,  some  temporaries  are  sparingly  used. 

2.3  The  HULL  Island  Concept 

The  BRL  version  of  HULL  has  the  capability  for  treating  rigid 
perfectly  reflecting  structures  by  allowing  the  user  to  specify  ensembles  of 
cells  in  the  mesh  that  are  themselves  perfectly  reflecting.  These  reflecting 
cells  are  called  islands. 

In  the  current  version  of  HULL,  a  cell  is  either  an  island  or  it  is 
not.  This  requires  that  all  structures  be  represented  as  combinations  of 
right  rectangular  solids.  Analyses  performed  by  the  AFWL  had  addressed  this 
approach,  specifically  for  the  MX  application  (reference  8).  They  found  that 
indiscriminate  zoning  would  tend  to  produce  too  high  an  overpressure  for 
certain  zones  adjacent  to  islands  when  the  islands  represented  a  slope  as  in 
Figure  la.  They  also  found  that  zoning  in  the  manner  shown  in  Figure  1b 
(i.e.,  where  they  required  the  center  of  the  air  cells  to  fall  on  the  hypo¬ 
thetical  surface  of  the  structure)  produced  good  agreement  with  theory  (in  the 
regular  reflection  region)  and  with  experimental  data  from  BRL. 

HULL  requires  information  on  each  of  the  six  faces  of  a  3-D  cell  at 
time  n  and  n+£.  It  derives  the  pressures  and  velocities  at  these  boundaries 
at  these  times  from  information  from  the  adjacent  cells  at  time  n. 

When  processing  a  given  cell  in  the  mesh,  the  left,  aft  and  bottom 
cells  have  already  been  updated  to  n+1;  the  information  is  saved  at  the 
boundaries  in  additional  arrays. 

Islands  use  reflective  boundary  conditions  to  represent  rigid  struc¬ 
tures.  A  reflective  boundary  is  one  where  the  normal  component  of  velocity  at 
the  surface  is  identically  zero.  Therefore  the  normal  component  of  velocity 
inside  an  island  cell  is  conceptually  equal  and  opposite  in  sign  to  the  normal 
velocity  component  at  the  same  distance  from  the  boundary  in  the  adjoining 
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a)  Early  Zoning  Method  for  a  Wedge 


b)  More  Recent  Zoning  Method 


Figure  1.  Zoning  a  Wedge 
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fluid  cell,  as  illustrated  in  Figure  2.  Using  linear  interpolation  and  noting 
that  the  velocity  on  the  boundary  joining  the  island  and  fluid  cell  should  be 
zero  yields: 


whence 


*(ui+1+u  )  =  0 


Section  3 
SHORE  CELLS 

This  section  describes  shore  cells  which  are  intended  for  use  in 
representing  rigid  structures  in  the  2-D  (Cartesian  and  cylindrical)  BRL  HULL 
code.  The  shore  cell  approach  allows  the  introduction  of  cells  into  HULL  that 
are  half  fluid  (normally  air)  and  half  perfectly  reflective.  Most  HULL  imple¬ 
mentations  currently  allow  cells  that:  (1)  consist  entirely  of  a  fluid  or  (2) 
are  perfect  reflectors  or  islands  (see  Section  2).  A  solid  structure  has  in 
the  past  been  simulated  by  an  ensemble  of  islands.  However,  not  all  struc¬ 
tures  can  be  well  represented  with  the  island  approach.  For  example,  a  wedge 
has  been  simulated  by  a  set  of  islands  arranged  in  a  stairstep  manner.  It  has 
been  shown  (see  Appendix  B  for  a  brief  review)  that  the  stairstep  simulation 
is  inadequate  for  certain  important  problems.  Therefore,  the  concept  of 
partially  reflective  cells  was  investigated  as  one  possible  solution  to  some 
of  the  problems  of  interest  to  the  BRL. 

In  order  to  improve  the  simulation  of  certain  structures  (e.g., 
ramps  or  other  surfaces  not  aligned  to  the  coordinate  system),  shore  cells 
have  been  implemented  in  2-D  BRL  HULL,  and  a  conceptual  approach  to  three 
dimensions  has  been  formulated.  A  shore  cell  in  two  dimensions  can  be 
imagined  as  a  cell  cut  in  half  along  its  diagonal.  (In  cylindrical  coordi¬ 
nates  the  'diagonal'  is  a  surface  of  revolution  which  cuts  the  cell  only 
approximately  in  half;  nevertheless,  in  the  ensuing  discussion  a  shore  cell  is 
referred  to  as  if  it  were  cut  in  half.)  Half  of  a  shore  cell  is  fluid  and  the 
other  half  is  perfectly  rigid  and  reflective,  or  in  HULL  terminology,  half¬ 
fluid  and  half-island. 

3.1  Introduction  to  Shore  Cells 


A  2-D  shore  cell  is  by  definition  a  cell  which  is  one-half  fluid 
and  one-half  island.  Furthermore,  the  boundary  between  these  two  halves  must 
be  one  of  the  two  cell  diagonals.  Therefore,  four  shore  cell  orientations  are 
possible  as  shown  in  Figure  3.  To  illustrate  how  the  diagonal  affects  the 
flow,  boundary  conditions  can  be  calculated  interior  to  the  shore  cell. 

Two  basic  assumptions  are  made  about  pressure  in  a  shore  cell:  (1)  the 
pressure  is  linear  parallel  to  the  diagonal  and  (2)  the  partial  derivative  of 
pressure,  on  the  diagonal,  normal  to  the  diagonal  is  zero.  These  lead  to  the 
computation  of  virtual  pressure  on  the  island  sides  of  shore  cells.  Consider 
Figure  4.  By  geometric  considerations 

LC  =  CL’  =  gj  =  -y.  sin  a  (16) 


B7?  =  CB  =  g2  =  ^  cos  a  (17) 
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3y  the  linearity  assumption 


PC  = 


(g2  PL  +  81  Pg) 


(g2  + 


i^T 


2  2 
cos  a  PL  +  sin  a  Pg 


(18) 


Where  P^,  P^,  and  P^,  are  the  pressures  at  points  L,  B,  and  C 

respectively.  Also  by  linearity, 


and 


or 


and 


pc  •  MPL  +  pl0 
PC  *  *  (PB*V> 


PL'  =  2  PC  -  PL 


P  •  =  2  P  =  P  . 

C  B 


(19) 

(20) 

(21) 

(22) 


By  reflectivity, 


and 


P  =  p  i  =  2  P  -  P 
R  PL  C  L 


P=P  :  2  P  -  P 
A  rB*  c  C  B 


(23) 

(24) 


Interpolated  results  are  only  expected  to  be  physically  meaningful  within 
the  confines  of  a  cell.  The  pressure  on  the  boundary  points  L  and  B  will  in 
the  physical  case  be  positive.  Experience  has  shown  that  when  properly  used, 
HULL  will  preserve  this  positivity.  However,  note,  that  while  interpolation 
to  point  L'  in  Figure  I  will  always  yield  a  positive  number,  it  is  possible 
that  the  pressure  at  B’  will  be  negative.  Point  B’  can  be  considered  to  be 
"outside"  of  the  cell,  and  its  value  is  only  a  numerical  convenience  to 
calculating  the  gradient  of  pressure  across  the  cell.  It  should  not  be  forced 
to  be  positive. 

It  is  also  assumed  that  the  component  of  velocity  normal  to  the  cell 
diagonal  must  vanish  at  the  cell  center.  This  is  imposed  at  the  end  of  each 
time  step.  Let 


-»■  +  •»• 

u  =  u  i  +  v  i 

be  the  computed  velocity  at  the  center  of  the  cell  where  and  "|j  are  the  usual 
unit  vectors  in  the  x  and  y  directions,  respectively.  One  can  also  write 


+  N  * 
u  =  u  n 


X 
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where  r?  and  t  are  unit  vectors  normal  to  and  tangent  to  the  diagonal, 
respectively.  By  the  assumptions,  for  the  orientation  in  Figure  4. 


uN  =  0 


and 

T 

u  =  -u  sin  a  +  v  cos  a 


then,  u  and  v  at  the  cell  center  are  recomputed  by 

u  =  -u?  sin  a 
v  =  uT  cos  a. 


(25) 

(26) 


(27) 

(28) 


3.2  Neighboring  Cells 


For  a  fluid  cell  calculation  without  shore  or  island  cells,  only  the 
fluid-fluid  interactions  need  to  be  considered.  Without  any  loss  in  gener¬ 
ality,  only  one  combination  of  two  cells  is  inspected.  This  approach  is  valid 
whether  sweeping  the  mesh  from  left  to  right  or  from  bottom  to  top  (or  from 
fore  to  aft  in  3-D).  When  islands  are  added,  four  possible  interactions  must 
be  inspected.  These  are: 


Interactions 


Boundary  between  cell  is 


1 .  fluid-fluid 

2.  fluid-island 

3.  island-fluid 

4.  island-island 


normal 

reflective 

reflective 

ignored 


The  order  chosen  to  calculate  the  cell  quantities  is  important  when 
handling  the  data.  Therefore,  four  interactions  were  listed  above  instead  of 
only  three.  The  additional  interactions  required  with  the  addition  of  shore 
cells  are: 


5. 

fluid-shore  (fluid) 

normal 

6. 

fluid-shore  (island) 

reflective 

7. 

shore  ( fluid) -fluid 

normal 

8. 

shore  (island) -fluid 

reflective 

9. 

island-shore  (fluid) 

reflective 

10. 

island-shore  (island) 

ignored 

11. 

shore  ( fluid )-island 

reflective 

12. 

shore  ( island)-island 

ignored 

13. 

shore  ( fluid) -shore  (fluid) 

normal 

14. 

shore  ( fluid )-shore  (island) 

reflective 

15. 

shore  ( island )-shore  (fluid) 

reflective 

16. 

shore  ( is land) -shore  (island) 

ignored 

a  shore  cell  can  have  its  fluid  side  or 

island  side  : 
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3.3  Shore  Cell  Modifications  to  the  Volume.  Density,  and 
Mass  Computations 

In  general,  the  most  important  change  for  shore  cells  is  in  the  correct 
expression  of  the  density.  The  usual  density  calculations  assume  a  full  cell. 
The  shore  cell  changes  generally  calculate  a  partial  cell's  volume  and  then 
subsequently  calculate  the  density  as  the  ratio  of  the  mass  of  fluid  in  the 
cell  to  its  volume.  This  approach  keeps  the  coding  from  obscuring  the 
physics. 

Numerous  minor  changes  are  required  for  the  density  computation.  Since 
these  changes  are  spread  all  over  the  code,  it  is  worthwhile  reviewing  their 
calculation  so  that  individuals  making  future  revisions  will  be  aware  of  what 
is  being  done. 

3.3.I  Modifications  of  Shore  Cell  Volumes.  Density  and  Mass  in 
2-D  Cartesian  Coordinates 


The  volume,  V,  of  a  2-D  Cartesian  cell  of  height,  Ay,  and  width,  Ax,  is 


and  mass  is 


v  =  Ax&y 


m  =  pV 


The  volume  of  fluid  in  a  2-D  Cartesian  shore  cell  is 


Vf  =  J  Ax  Ay  =  JV 

and  mass  is 


m  =  i  p  V 


3-3.2  Modification  of  Shore  Cell  Volume.  Density  and  Mass  in 
2-D  Cylindrical  Coordinates 

In  cylindrical  coordinates  the  volume,  V,  of  a  cell  of  height,  Ay,  and 
width,  Ax  is  that  for  a  toroidal  cell 


V 


Ax 

2 


1  Ax]  2 

c  "  2 


it  Ay 


or  V  =  2>r  RQ  Ax  Ay 
where  x  is  the  radial  component. 
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For  cylindrical  coordinates,  HULL  defines  PI  =  tt  ,  RC^  =  J(x^  +  xi+i) 

and  TAU^  =  J  PI  RC*  Axj..  An  obscurity  exists  in  HULL  to  allow  use  of  the  same 
code  for  Cartesian  geometry  as  for  cylindrical  geometry.  The  variable  PI  is 
set  to  £  and  Rc  to  1  so  that  in  Cartesian  geometry  TAU  =  Ax.  In  either 
geometry 


V  =  TAU  Ay. 


No  attempt  has  been  made  to  remove  this  obsurity,  and  furthermore, 
such  use  of  TAU  may  have  been  included  in  coding  for  SHORE  cells.  Eventually, 
such  coding  should  be  removed  from  HULL  for  although  it  can  be  considered  to 
be  computationally  efficient  in  some  sense,  its  use  could  easily  lead  to 
future  errors. 

Consider  the  shore  cells  shown  in  Figure  5.  The  number,  L  contained 
within  the  fluid  portion  of  the  shore  cell  identifies  its  orientation.  If  L 
is  greater  than  2  the  volume  is  greater  than  one-half  the  volume  of  the  cell; 
otherwise  it  is  less. 

In  general,  the  volume  of  a  solid  of  revolution  can  be  calculated  as 
shown  in  Figure  6.  From  the  Theorem  of  Pappus  the  volume  of  revolution  of  a 
cell  can  be  obtained  easily  if  the  center  of  gravity  and  cross-sectional  area 
are  known.  The  center  of  gravity  on  the  fluid  side  is  one-third  Ax  and  one- 
third  Ay  from  the  boundaries  adjacent  to  the  fluid.  Therefore,  the  volume  of 
fluid  is 


Vf  =  tt(Rq  ±  ^  Ax)  Ax  Ay 


instead  of  2irRc  AxAy.  The  sign  depends  on  which  side  of  the  diagonal  is  being 

considered.  In  other  words,  the  radius  to  the  centroid  is  no  longer  Rc  but 

instead  is  Rc  ±  Jl  Ax  and  the  cross-sectional  area  has  been  cut  in  half.  In 
6 

our  implementation  an  array  TAUS  has  been  defined  that  stores  the  difference 
that  needs  to  be  added  to  (or  subtracted  from)  TAU.  Its  value  is 


TAUS 


T  <Ax)‘ 


(29) 


so  that  in  cylindrical  coordinates  the  volume  of  a  shore  cell  of  type  L>2  is 


Vf  =  (  jTAU+TAUS)Ay 
For  type  L  <  2, 
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Figure  5. 

Shore  Cell  Nomenclature. 
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Figure 

6.  Solid  of  revolution 

Vf  =  (J  TAU-TAUS )A y 


3.H  Shore  Cell  Modifications  to  the  Lagrangian  Phase 

This  section  provides  an  overview  of  the  physics  of  the  implementation  of 
the  shore  cell  concept  for  the  Lagrangian  Phase.  The  approach  does  not 
require  extensive  modification  to  the  code's  architecture,  and  in  fact,  has 
been  implemented  with  an  appropriate  POST  option  so  that  unless  the  user 
explicitly  invokes  the  shore  cell  option,  its  source  code  will  not  be  included 
in  HULL. 


The  conditions  on  fluid  boundaries  of  a  shore  cell  are  computed  in  the 
normal  manner.  That  is,  if  a  cell  adjacent  to  the  fluid  boundary  is  a  fluid 
cell  or  the  fluid  side  of  a  shore  cell,  then  the  boundary  is  treated  the  same 
as  two  adjoining  fluid  cells.  If  the  neighbor  is  an  island  or  the  reflective 
side  of  a  shore  cell,  then  the  boundary  is  treated  as  reflective. 


The  difference  equation  analog  of  Equation  (2)  used  in  HULL  for  2-D 
Cartesian  fluid-cells  is 


n+1  n  .  .  (Pi-f,j  ”  Pi+f,j)  Ayj 

u.  .  =  u.  .  +  At  - - “■ 


*  .  «  —  v* .  , 

l,j  i.j 


n 
m.  . 
i»3 


and 


„  .  (XL  -  pZlib’ 

V.  .  =  V.  .  +  At  — 2 - 


1,J  i,J 


n 

m.  . 


-  At  G 


(30) 


(3D 


J 


Here  Gj  is  the  "gravity  potential"  which  prevents  any  change  in  Vjj  in  an 
undisturbed  ambient  atmosphere.  Note  that 


Ay 

n 


m 


ij 


Ax.^  Ay^ 

Ax  nr  . 
i  iJ 


Ax.  V.  .  p‘‘ 
i  ij  ij 


1 


Ax.  pn 
i  ij 


(32) 


and  similarly 


1 


ij 


IJ 


Ay.  Pn, 
i  ij 


(33) 


where  V^j  is  the  volume  of  the  cell. 
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By  a  "clever"  choice  of  multipliers  the  same  coding  is  used  for  both 
Cartesian  and  cylindrical  geometry.  In  either  case  1/p|?j  is  represented  by 


For  shore  cells  the  volume  of  fluid  and  density  in  the  cell  are 
explicitly  computed.  The  pressures  on  the  two  fluid  sides  of  the  shore  cell 
at  time  n+J  are  computed  as  for  cells  that  are  entirely  fluid.  Assuming  the 
orientation  described  in  Figure  4,  Pl  and  Pg  are  computed  in  the  normal  manner 
and  Pr  and  P^  are  computed  as  described  in  Section  3.1. 


Then 


and 


n+1 

'i, 


+  At 


n+1 

ViJ 


+  At 


4yj  pi,j 


4t  Gj 


(34) 


(35) 


Since  these  equations  do  not  depend  on  mass  or  volume,  they  are  applica¬ 
ble  to  both  Cartesian  and  cylindrical  geometry  and  to  both  all  fluid  cells  and 
shore  cells. 

The  energy  equation  (Equation  3)  in  Lagrangian  coordinates  can  be 
rewritten  as: 

+  IV  •  (pu)  +V6  (36) 

at  p 


When  formulating  a  difference  equation  involving  the  divergence  it  is  conven¬ 
ient  to  write  down  its  definition  and  then  work  from  there.  The  divergence 
used  in  the  HULL  difference  equations  is  the  "average"  value  in  the  cell.  For 
a  vector  "F, 


(average  v.t)  f  dV  =  f„  v"fr  dV, 

J  V  J  V 

By  Gauss'  theorem 

JvvtdV  =  Jst.dt 

where  V  is  volume  of  the  cell,  S  is  the  surface  area  of  the  cell,  and  d"X  is 
the  product  of  the  differential  of  area  with  the  unit  normal  on  the  surface. 


"or  a  rectangular  cell  of  size  Ax  and  Ay  (with  Az  =  1),  the  divergence 

~y 

of  pu  is  approximately 

.  -Ay(pu).  +  Ay(pu)_  -  Ax(pv)  +  Ax(pv) 

V  •  (pu)  =  - ± L 2 - *  (38) 

where  (pu)^  means  its  value  on  the  left  side  of  the  cell,  R  =  right,  B  = 
bottom,  and  A  =  top.  Remember,  u  is  the  velocity's  component  along  x  and  v  is 
its  component  along  y. 


Substituting  in  Equation  (36)  with  m  =  pV. 

AE  =  ^  {|(pu)L  -  (pu)R]  Ay  +  |(pv)B  -  (pvjjjAx  -AtvG 


The  same  sort  of  argument  for  cylindrical  coordinates  produces 

AE  =  jjp-  ||(xpu)L  -  (xpu)R|  27tAy  +  [ (pv)Q  -  (pv)A|  211  Rc  Ax  j 

-  At  vG. 


(39b) 


A  shore  cell  on  the  other  hand  has  three  sides  enclosing  the  fluid 
instead  of  four  (in  two-dimensions).  The  contribution  to  the  divergence 
across  the  diagonal  is  zero  and  therefore  only  the  sides  of  the  cell  contri¬ 
bute.  Furthermore,  since  the  velocities  at  the  reflective  sides  of  the  shore 
cell  are  zero,  the  original  energy  equation  could  be  used  if  it  were  not  for 
the  gravity  potential  term.  The  gravity  potential  term  computed  in  HULL  for  a 
full  cell  is 


where  the  pressures  and  densities  are  taken  from  an  ambient  atmosphere  and 
computed  just  as  they  would  be  for  the  difference  equations.  For  all  fluid 
cells,  for  both  Cartesian  and  cylindrical  geometry,  the  same  term  prevents  a 
change  of  vertical  velocity  and  accounts  for  the  effect  of  gravity  on  energy 
in  an  ambient  atmosphere.  The  Gj  stored  is  correct  for  the  momentum  equation 
for  shore  cells  in  both  cylindrical  and  Cartesian  coordinates.  However,  we 
must  multiply  the  gravity  term  Gj  by  V/V^  in  the  energy  equation.  For 
Cartesian  coordinate  shore  cells 

ae  =  ^  |[(pu)l“(pu)r]  +  f( Pv Pv ^ A J  Axj>-  2  At  vG 
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For  cylindrical  shore  cells 


AE  =  ^  {[<xpu)L-  (xpu)^27rAy 

+  ^{(pv)B-(pv)j2TrRcAx| 

2u  R  AxAy 

- n -  At  vG 

f 

Here,  2ttRc  AxAy  =  Ay  TAU  is  the  volume  of  the  cell  and  Vf  is  the  volume  of  the 
fluid  part  of  the  cell.  This  will  have  been  computed  to  find  density  from 
mass. 


3-5  Shore  Cell  Modifications  to  The  Fluxer  Phase 

The  approach  taken  for  fluxing  mass,  momentum,  and  energy  is 
similar  to  the  original  HULL  approach.  However  if  the  donating  cell  is  a 
shore  cell,  the  density,  p,  is  computed  and  the  mass  flux  computed  using  p 

explicitly, 

Mass  flux  =  P(uAt)  A 


where,  u  is  the  velocity  component  perpendicular  to  the  cell  boundary,  A  is 
the  cell  boundary  area,  and  At  is  the  time  step. 

After  fluxing  mass,  momentum  and  total  energy,  the  boundary 
conditions  are  applied  so  the  velocity  normal  to  the  reflecting  surface  will 
be  zero  at  the  cell  center  (see  Section  3.1).  This  is  done  while  strictly 
conserving  mass  and  total  energy.  The  result  when  the  flow  does  not  happen  to 
come  out  strictly  parallel  to  the  diagonal  wall  is  to  convert  kinetic  energy 
into  internal  energy. 


Section  4 

ADDITIONAL  DETAILS  OF  THE  IMPLEMENTATION  OF  SHORE  CELLS 

4.1  Modifications  to  PLANK 

The  changes  to  PLANK  primarily  involved  the  modifications  necessary  to 
include  an  additional  option,  SHORE,  into  the  HULL  system.  At  the  same  time 
the  OPT  array  was  reduced  to  what  it  had  been  when  HULL  was  initially 
installed  at  BRL. 

4.2  Modifications  to  KEEL 


Appropriate  comments  were  added  to  KEEL,  the  grid  generator  for  HULL,  to 
maintain  the  current  state  of  documentation;  these  should  be  self  explanatory. 
The  option  SHORE  was  also  included  in  the  HULL  z-block,  since  each  program  in 
HULL  should  know  of  its  existence,  and  it  would  be  inappropriate  to  exclude  it 
from  the  restart  files. 

The  architecture  of  KEEL  has  not  evolved  to  an  easily  modifiable  or 
maintainable  state,  and  it  was  out  of  the  scope  of  this  effort  to  change  this 
circumstance.  Therefore,  a  technique  was  developed  that  avoided  a  significant 
rewrite;  unfortunately,  the  approach  is  obscure  to  the  casual  user.  For  that 
reason,  some  detailed  comments  are  appropriate  even  though  they  will  not  make 
KEEL  easy  to  understand. 

Essentially  KEEL  is  designed  for  convenience  of  use;  it  allows  simple 
descriptions  of  geometric  objects  or  regions  to  be  specified  by  the  analyst 
and  then  will  assign  hydrodynamic  values  to  each  cell  by  mass-weighting  the 
hydrodynamic  values  for  subcells  (discussed  below).  The  allowable  two-dimen¬ 
sional  objects  or  regions  are:  rectangles,  triangles  and  circles.  The 
allowable  three-dimensional  objects  or  regions  are:  boxes,  tetrahedrons, 
spheres,  cylinders  and  cones.  The  standard  2-D  HULL  permits  multi-material 
input.  The  3-D  HULL  (or  the  2-D  shore  option)  restricts  the  contents  of  any 
region  to  be  air  (the  state  of  the  air  may  be  different  in  each  region)  or  an 
unyielding  solid,  designated  as  ISLAND. 

Let  it  suffice  to  say  that  very  complex  geometrical  shapes  can  be  speci¬ 
fied  rather  easily,  since  objects  can  be  added  or  deleted  as  needed.  The 
coding  to  perform  this  in  KEEL,  however,  isn't  straightforward.  This  is 
partially  the  result  of  many  years  of  disjoint  development. 

All  the  regions  containing  fluid  are  processed  before  the  ISLAND  regions. 
Currently,  for  the  fluid  regions,  each  cell  is  divided  into  subcells  with  3 
partitions  in  each  direction.  This  produces  9  subcells  for  2-D  geometry  and 
would  produce  27  subcells  for  3-D.  If  the  center  of  a  subcell  is  in  a  region, 
that  subcell's  portion  of  the  cell  is  assigned  the  hydrodynamic  values  of  that 
region. 

The  island  regions  are  processed  last.  Each  potential  island  cell  is 
considered  as  a  unit  (one  subcell).  If  the  center  of  the  cell  is  inside  anv 
island  region  the  cell  is  designated  an  island  cell. 

A  final  check  fills  the  unfilled  portions  of  cells  with  ambient  air. 


If  the  center  of  a  subcell  should  lie  on  the  common  boundary  between  two 
fluid  regions,  that  subcell  might  be  included  in  both  regions.  In  this  case 
the  cell  might  be  overfilled  and  the  program  would  halt  with  a  warning  print. 
Conversely,  it  might  be  excluded  from  both  regions.  In  this  case  the  omitted 
volume  would  be  filled  with  ambient  air.  Overlap  of  island  regions  with 
either  fluid  or  other  island  regions  does  not  cause  trouble  because  the 
islands  are  entered  last  and  entire  cells  are  either  island  or  not. 

For  the  2-D  shore  cell  option  the  fluid  regions  were  treated  just  as 
before  except  there  are  4  partitions  in  each  direction.  This  gives  16 
subcells.  This  partition  into  16  subcells  is  retained  for  the  rigid  regions 
designated  by  the  word  SHORE  (or  ISLAND).  If  any  subcell  in  the  cell  is  in  a 
shore  region,  the  appropriate  bit  in  a  16  bit  true-false  piece  of  a  word  is 
set  true.  After  all  the  regions  have  been  processed,  this  true-false  informa¬ 
tion  is  used  to  determine  whether  the  cell  should  be  fluid,  island,  or  one  of 
the  four  possible  shore  cells. 

This  determination  of  the  type  of  cell  is  carried  out  as  the  last  part  of 
the  final  pass  that  fills  the  unfilled  portion  of  cells  with  ambient  air.  The 
coding  includes  adjustments  to  the  content  of  cells  to  the  proper  level  for 
shore  cells  or  for  discarded  island  subcells. 

If  there  are  less  than  six  subcells  filled  with  island  material,  then  the 
cell  is  considered  to  be  entirely  fluid  and  the  subeells  that  were  islands  are 
filled  with  fluid.  If  there  are  more  than  thirteen  subcells  that  are  filled 
with  island,  then  the  entire  cell  is  made  an  island. 

The  remainder  of  the  cells  are  treated  as  potential  shore  cells.  Depend¬ 
ing  on  the  outcome  of  the  following  tests  the  cell  may  become  a  shore  cell,  an 
island,  or  a  fluid  cell.  These  potential  shore  cells  are  checked  to  see  if 
their  "corner"  subcells  are  island  or  not.  If  the  six  corner  subcells  marked 
with  an  X  in  Figure  7a  contain  island  material  then  the  cell  can  become  a 
shore  cell  with  an  orientation  tag  of  1 .  This  is  recognized  by  setting  a 
logical  flag,  LSI,  to  TRUE.  A  similar  test  is  performed  for  each  of  the  other 
orientations  and  corresponding  logical  flags  are  set.  There  are  three  possi¬ 
bilities: 

1.  If  LSI,  LS2,  LS3,  and  LS4  are  all  false,  the  cell  Is  fluid. 

2.  If  any  two  or  more  of  LSI,  LS2,  LS3,  and  LS4  are  true,  the  cell  is 
an  island. 

3.  If  exactly  one  of  LSI,  LS2,  LS3,  and  LS4  is  true,  say  LSi,  then  the 
cell  is  a  shore  cell  of  type  i. 

Note  that  the  cell  in  Figure  4b  does  not  pass  the  corner  criterion  for 
any  shore  orientation,  and  therefore  becomes  a  fluid  cell.  It  is  unlikely 
that  this  case  will  be  encountered  except  possibly  at  point  junctions  of  three 
structures,  in  which  case  making  it  a  fluid  could  be  appropriate. 

Some  simplifying  assumptions  have  been  made  to  avoid  most  of  the  changes 
that  would  be  necessary  at  the  mesh  boundaries.  For  the  time  being,  shore 
cells  are  not  allowed  to  be  adjacent  to  a  boundary  except  at  either  a  reflec¬ 
tive  left  boundary  (LBOUND  =0),  or  a  reflective  bottom  boundary  (BB0UND=0). 
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The  remaining  changes  to  KEEL  involve  calculating  the  volumes  of  the 
cells  and  subcells  and  modifying  the  materials  maps. 

4.3  Modifications  to  HULL:  Adding  the  SHORE  Option 

The  modifications  to  HULL  in  IDENT  HULBS1  add  the  shore  option  to  the 
code  and  properly  treat  its  value  in  the  z-block.  Also,  comments  were  added 
consistent  with  the  practice  at  BRL. 

4.4  Modifications  to  HULL:  2-D  Consideration  and  Lagrangian  Phase 
Modifications 


Information  about  a  shore  cell  is  kept  in  the  H-array  as  another  hydro- 
dynamic  variable.  This  is  consistent  with  Vector  HULL  architecture,  and  is 
the  appropriate  way  to  introduce  shore  cells  to  make  it  possible  to  vectorize 
the  code. 

This  additional  hydrodynamic  variable  contains  the  shore  orientation 
indicated  by  a  real  value  of  1.0,  2.0,  3-0,  or  4.0.  If  this  variable  is  0.0 
the  cell  is  fluid,  and  if  -1.0  an  island.  Rather  than  use  a  real  variable  for 
comparison,  an  integer  LSX  is  used  where  X  can  be  I  for  the  current  cell,  R 
for  the  right  cell,  or  A  for  the  above  cell — consistent  with  HULL  mnemonics. 

A  change  made  to  the  time  step  calculation  was  included  in  this  update. 

In  HI  (the  Lagrangian  phase)  it  was  also  necessary  to  consider  how 
density  should  be  computed  at  a  boundary  between  two  cells,  one  or  both  of 
which  are  shore  cells. 

As  an  example  of  the  approach  used,  consider  a  fully  fluid  cell  with  an 
LSA=4  shore  cell  above  it  (its  fluid  side  faces  the  fully  fluid  cell). 
Consider  the  case  for  cylindrical  geometry. 

In  HULL  the  mass  of  two  all-fluid  cells  is 
AMA  =  H(N1+5 )  +  H(NA1+5) . 

The  density  at  that  boundary  is  calculated  as 

RHOA  =  AMA/(TAU(I)  »  (DY(J)  +  DY(J+1))) 

In  order  to  calculate  the  density  with  this  equation  we  choose  to  modify  the 
mass  associated  with  each  shore  cell  at  the  boundary. 

If  Mf  represents  the  mass  of  fluid  in  the  shore  cell,  and  Vf  the  volume 
of  fluid,  and  m  represents  the  mass  of  the  cell  if  filled  with  the  fluid  (the 
quantity  desired),  and  V  its  volume,  we  can  write: 


or 

^f _ 

“  '  (Vf/V) 
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Then,  for  a  shore  cell  above  a  fluid  cell,  the  apparent  mass  for  this  density 
interpolation  is 


AMA  =  H(N1+5)  mf/(Vf/V)  . 

For  the  chosen  case,  i.e.,  LSA  =  4, 

V  =  TAU(I)  «  DY( J+1 ) 


and 


or 


Vf  =  (Jj  TAU ( I )  +  TAUS(I))  *  DY(J-fl) 

Vf  .  TAUS(I) 

V~  ~  *  +  TAU(I) 


Similarly  for  LSA  =  1 

If  ,  TAUS(I) 
V  =  *  ‘  TAU(I) 


Section  5 

PRELIMINARY  RESULTS 


5.1  Initial  Code  Checking 

In  late  May  1981,  SAI's  files  containing  the  modifications  for  shore 
cells  were  transferred  to  the  BRL's  CDC  7600  computer  and  code  checking  was 
begun.  The  majority  of  the  test  computations  were  performed  by  the  BRL.  A 
number  of  changes  were  made  and  tested.  A  revised  set  of  shore  modifications 
files,  with  additional  improvements,  was  sent  from  SAI  to  BRL  in  mid- 
September  . 

Much  of  the  checking  of  the  shore  changes  was  accomplished  by  review  of 
the  coding  and  considering  the  correctness  under  all  permissable  conditions. 
This  is  not  completely  reliable,  but  is  more  feasible  than  actually  testing 
all  possible  code  configurations. 

5.2  Test  Computations  (2-D  Cartesian  Coordinates) 

Testing  is  still  continuing  and  will  be  documented  in  a  later  BRL 
report.  To  illustrate  the  use  and  effect  of  shore  cells,  three  computations 
from  the  first  series  of  tests  will  be  briefly  discussed  here.  All  three  used 
a  100  by  100  grid  of  square  cells  5  cm  on  a  side,  had  the  same  ambient  atmos¬ 
phere,  and  had  the  same  68.95  kPa  (10  psi)  step  shock  striking  a  square 
unyielding  target  at  an  angle  of  45°  to  the  front  sides.  The  runs  were 
initiated  with  the  shock  front  20  cm  from  the  leading  corner  of  the  square 
target.  The  initial  time  was  set  at  4.1  ms  and  the  computations  ran  to  16  ms. 

Problem  181.0714,  which  we  call  the  shore  diamond  test,  used  the  shore 
modifications.  The  step  shock  was  input  from  the  left  boundary,  and  struck 
the  square  target  (the  diamond)  at  45°  to  the  front  edges  (See  Figure  8). 
Problem  181.0715,  which  we  call  the  oblique  shock  test,  had  none  of  the  shore 
modifications.  The  step  shock  moved  from  the  lower  left  toward  the  upper 
right  at  45°  to  the  square  computational  field  and  the  square  target,  as  indi¬ 
cated  in  Figure  9.  Problem  181.1028,  which  we  call  the  staircase  diamond 
test,  was  a  duplicate  of  the  shore  cell  diamond  test  except  the  shore  changes 
were  not  used.  The  target  was  compressed  slightly  so  the  shore  cells  in 
Problem  181.0714  were  replaced  by  all  fluid  cells. 

We  include  two  of  the  HULL-produced  plots  from  each  of  these  runs. 
Figures  8,  9,  and  10  are  plots  of  the  vector  velocities  at  14  ms  for  the  shore 
diamond  test,  the  oblique  shock  test,  and  the  stairstep  diamond  test, 
respectively.  This  time  was  selected  to  point  out  the  strong  vortices  that 
form  behind  the  targets.  A  close  inspection  of  Figures  8  or  10  reveals  an 
apparent  asymmetry  in  the  velocity  vectors.  This  apparent  asymmetry  in  the 
vicinity  of  the  target  is  entirely  due  to  the  plotting  procedure  which  plotted 
velocity  vectors  for  cells  in  alternating  rows  and  columns.  The  asymmetry 
near  the  bottom  boundary  relative  to  positions  near  the  top  boundary  is  real. 
BRL  HULL  does  not  have  a  satisfactory  transmissive  lower  boundary.  Since  we 
wanted  to  compare  results  with  the  oblique  shock  test,  which  is  transmissive 
parallel  to  the  shock  front,  we  made  the  lower  boundary  reflective  and  the 
upper  boundary  transmissive.  The  effect  of  this  near  the  target  was  negligi¬ 
ble  at  14  ms.  The  plotting  from  alternate  cells  in  the  oblique  shock  test  is 
symmetric. 
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Figures  11,  12,  and  13,  are  overpressure  contour  plots  at  10  ms  for  the 
same  three  HULL  problems  (in  the  same  order  as  before).  This  particular  time 
was  chosen  because  the  plotting  procedure  selected  the  same  contour  levels  for 
all  three  cases.  The  apparent  outlines  of  the  target  are  an  artifice  of  the 
plotting  program  (A  very  low  pressure  is  assigned  to  the  target  cells  and  the 
interpolation  produces  several  close  contours .)  There  has  been  no  attempt  to 
modify  the  plotting  procedure  for  shore  cells.  The  pressure  contour  plotting 
procedure  assumes  square  cells,  even  to  the  point  of  recomputing  the  pressure 
in  a  cell  assuming  the  entire  cell  is  fluid.  This  does  no  real  harm  as  long 
as  the  observer  remembers  that  the  apparent  target  outline  is  severely 
distorted. 

A  definitive  comparison  of  results  from  these  runs  is  not  possible  from 
these  plots.  Our  intention  was  to  compare  values  at  various  "stations"  around 
the  perimeter  of  the  target.  The  values  at  a  station  are  the  values  for  the 
cell  in  which  the  station  is  located.  The  results  for  the  shore  diamond  test 
and  the  oblique  shock  test  were  not  quite  the  same.  This  is  largely  because 
the  cells  and  the  centers  of  the  cells  were  oriented  differently  with  respect 
to  the  target.  The  centers  of  the  shore  cells  are  on  the  target,  while  the 
centers  of  bordering  cells  for  the  other  case  were  2.5  cm  away  from  the  target 
edge.  An  examination  of  overpressure  plots  (not  shown  here)  at  various  points 
around  the  target  indicates  very  similar  results  on  the  front  (windward) 
sides,  quite  different  results  just  around  the  corners,  and  fair  agreement 
further  along  the  back  side. 

Figure  14  shows  the  history  of  average  overpressure  in  cells  whose 
sides,  or  diagonals,  form  the  front  of  the  target  and  average  overpressure  in 
cells  whose  sides,  or  diagonals,  form  the  back.  The  higher  curves  are,  of 
course,  from  the  front.  Occasional  squares  mark  the  shore  diamond  test 
results,  triangles  identify  the  off  angle  results,  and  the  stairstep  diamond 
results  have  no  superposed  symbol.  The  results  for  the  shore  diamond  test  and 
the  oblique  shock  test  are  close  enough  that  their  differences  may  be  due  to 
the  relative  locations  of  the  cell  centers.  The  average  pressure  for  the 
staircase  diamond  run  is  significantly  different.  Note  particularly  that  the 
difference  in  pressure  on  the  front  and  the  back  is  much  smaller  after  about 
10  ms. 


If  results  from  the  oblique  shock  test  are  accepted  as  correct  (this 
option  has  been  well  checked),  a  slanting  wall  of  shore  cells  is  significantly 
better  than  a  stairstep  wall  for  predicting  loads. 

5.3  Test  Computation  (2-D  Cylindrical  Coordinates) 

A  sample  run  with  cylindrical  coordinates  compiled  and  ran  with  no 
obvious  errors.  The  run  modeled  a  step  shock  moving  down  a  cylindrical  tube 
with  a  constricting  section  and  an  expanding  section.  Although  the  results 
appear  reasonable,  the  model  is  artificial  and  hence  this  case  does  not  serve 
as  a  check  for  correctness.  (In  fact,  a  minor  error  in  the  cylindrical 
coordinate  coding  was  later  found  and  corrected). 

A  cylindrical  shock  tube  test  is  planned  at  the  BRL  in  the  near  future. 
Data  from  this  test  will  be  used  to  check  the  shore  coding  for  cylindrical 
coordinates. 
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Figure  12.  Pressure  Contours  for  Oblique  Shock  Case 
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Figure  14.  Comparisons  of  average  pressures. 


Section  6 

CONCLUSIONS  AND  RECOMMENDATION 

The  results  of  the  tests  briefly  discussed  in  the  previous  section,  and 
other  preliminary  results,  suggest  that  the  shore  cells  will  prove  to  be  a 
useful  addition  to  HULL.  The  use  of  shore  cells  should  increase  the  number  of 
shapes  that  can  be  satisfactorily  modeled. 

It  should  be  pointed  out  that  not  all  shapes  can  be  modeled.  The  edges 
of  the  modeled  solid  will  follow  cell  boundaries  or  cell  diagonals.  A  very 
shallow  or  very  steep  ramp  cannot  be  modeled  without  introducing  cells  with 
very  high  aspect  ratios.  Use  of  such  cells  is  generally  unacceptable. 
Further,  a  sudden  change  of  slope,  say  from  30°  to  *45°,  would  also  cause 
modeling  difficulties.  To  maintain  smoothness  a  sudden  change  in  all  dimen¬ 
sions  would  be  necessary,  which  can  lead  to  other  inaccuracies.  Nevertheless, 
for  these  cases  a  more  realistic  model  can  be  formed  by  using  some  shore  cells 
than  without  them. 

Although  we  are  pleased  with  the  shore  cell  coding  for  2-D  Cartesian 
geometry,  some  further  testing  for  simplification  or  improvement  is  in  order. 
For  example,  it  is  not  definite  that  stagnation  (forcing  the  velocity  in  shore 
cells  to  be  parallel  to  the  diagonal)  is  needed.  If  one  considers  the  hydro- 
dynamic  variables  to  be  values  at  the  cell  center,  there  should  not  be  any 
velocity  component  normal  to  the  diagonal  in  shore  cells.  However,  if  the 
hydrodynamic  variables  are  values  associated  with  the  cell  center  (e.g., 
average  values  in  the  cell),  a  small  velocity  component  normal  to  the  cell 
diagonal  may  be  acceptable.  A  cursory  study  of  using  stagnation  versus  not 
using  it  was  effected  for  a  step  shock  striking  a  ramp.  The  only  significant 
result  apparently  was  to  produce  slightly  higher  peak  pressures  when  stagna¬ 
tion  was  not  used. 

At  the  time  of  this  writing,  the  accuracy  of  results  from  shore  cells 
with  2-D  cylindrical  coordinates  has  not  been  verified. 

In  view  of  the  apparent  success  of  2-D  shore  cells,  we  recommend 
proceeding  with  implementation  of  3-D  shore  cells  in  BRL  HULL. 
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Appendix  A 


DERIVATION  OF  TIME  DERIVATIVE  OF  PRESSURE 


This  appendix  presents  the  derivation  of  the  time  derivative  of 
the  pressure  in  order  to  demonstate  the  validity  of  the  differential 
equations  that  are  differenced  in  the  HULL  code. 


Substituting  equations  (4)  from  the  main  text  into  (3)  and  using 
(2)  one  can  show  that 


+  p  (V»u)  =  0 


(Al) 


and  then  if  one  defines  Yeff.  =  1  +  jj^r  it  also  can  be  shown  (using  (1))  that 

(A2) 


3t  *  <>  (7'“)  *  0 


The  proofs  follow.  Substituting  equation  (4)  into  (3)  yields 


P^T  (I  +  iu*u)  +  V*  (pu)  =  fju*g 


Now  since 

V*(pu)  =  (Vp)«u  +  p(V«u) 
it  follows  that 

P^T  (I  +  Ju*u)  +  (Vp)»u  +  p(V*u)  =  pu-g 
or  P^  +  jp^(u-u)  +  (Vp)*u  +  p(V*u)  =  p  a*  g 
or  p^  +  p  (V*u)  +  u  •  {p^  +  Vp  _  p^}  =  o 
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from  equation  (2)  in  the  main  text  it  follows  that 


+  p(7*u)  =  0 


which  is  the  first  thing  to  be  shown. 


Defining  Ygff  =  1  +  -^  and  writing  it  simply  as  Y 


3?  * 


1  dp  2  «  ]  dp 

*’  p  dt  "  I  dt  +  P  dt 


or  £1 

or  p  dt  *at  dt 


dl 


Substituting  the  above  for  from  equation  (Al)  yields 


dp 


Substituting  into  this  for  ^  f roo  (1): 


—  ^  +  IP(Vu)  +  p(  V«u)  =  0 

p  dt 


+  (p+IP)(V*u)  s  0 


or  ^  +  pY(V*u)  =  0 
dt 


which  was  the  second  thing  to  be  shown 


APPENDIX  B 


THE  CHANGE  DECK  TO  IMPLEMENT  SHORE  CELLS  IN  TWO-DIMENSIONAL  HULL 


Appendix  B 


This  appendix  contains  a  listing  of  the  change  deck  that  is  sent  to  the 
CDC  UPDATE  facility.  These  changes  represent  the  recommended  changes  to  the 
BRL  HULL  system  to  effect  shore  cells  and  were  developed  during  this  effort. 
In  addition,  various  changes  were  introduced  to  correct  errors  unrelated  to 
shore  cells  that  were  found  by  chance  during  the  implementation  phase  of  this 
contract.  These  recommended  changes  are  believed  to  be  correct,  however, 
additional  work  needs  to  be  done  to  fully  check  out  their  correctness. 

As  a  note  of  interest,  the  first  author  would  like  to  point  out  that 
the  typed  listing  in  this  appendix  is  produced  directly  from  the  magnetic 
media  where  the  changes  were  stored  while  being  developed.  With  current  tech¬ 
nology  we  can  produce  listings  of  letter  quality  without  the  usual  introduc¬ 
tion  of  errors  associated  with  retyping.  Furthermore,  the  changes  when  ready 
to  be  tested  were  sent  telephonically  to  BRL,  sent  back  to  the  sending  site, 
and  compared  at  the  sending  site,  character  by  character,  to  mitigate  communi¬ 
cation  errors. 

Finally,  the  ability  to  develop  the  changes  off-line  greatly  reduced 
associated  communications  costs. 


HULL  CHANGE  DECK  FOR  BRL:  SAI  SHORE  ISLANDS.  B.  CHAMBER 
<<<  31-OCT-81  >>>  VERSION  FOR  PLANK,  PROLOGUE,  AND  KEEL 
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SKIPPING  SETTING  OF  ISLANDS 
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ISHOR=0 
VOLTOT=0 . 0 

♦COPY  KEEL, KEEL. 1818, KEEL. 1821 
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♦LABEL  ENDSHORE 


VOLR=VOLR-TAUS (1+1) *DY (J) 

SINCE  THIS  IS  THE  LEFT-HAND  SIDE  OF  A  SHORE (AIR) 
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CELL  TO  THE  RIGHT  IS  SHORE (AIR)  OR  AIR 
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*COPY  HULL, HULL. 4052, HULL. 4070 
GO  TO  2000 
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2000  CONTINUE 

IF(H(NAl+4) .LE.O)  GO  TO  2500 

IF( (LSA.EQ.2) .OR. (LSA.EQ.3) )  GO  TO  2500 


CELL  ABOVE  IS  SHORE (AIR)  OR  AIR 
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* LABEL  ENDG2 

RHOA=H (Nl+5) /VOL 
♦COPY  HULL, HULL. 4156, HULL. 4158 
3000  CONTINUE 

IF (ABS (UR) .LE.VMIN)  UR=0 


EFFORT,  AND  SINCE  IT  HAS  NOT  YET  BEEN  TESTED,  IT  HAS 
NOT  BEEN  INCLUDED  IN  THIS  LISTING  OF  THE  CHANGES. 
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♦INSERT  HULL. 9226 
♦LABEL  ENDG2 

*/  HULL  CHANGE  DECK  FOR  BRL:  SAI  SHORE  ISLANDS.  B.  CHAMBERS 
*/  <<<  31-OCT-81  >>>  VERSION  FOR  HULL  PART 2  &  STATIONS 
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NOTE  ADDED  AFTER  REVIEW  OF  FINAL  DRAFT: 

THE  CURRENT  (FEB  82)  BRL  VERSION  HAS  BEEN  REVISED  TO 
"END  NOTE".  SINCE  THE  CHANGE  WAS  NOT  PART  OF  THIS 
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MULTIPLYING  BY  DY  IS  INTENTIONAL  BSC 3 ,  ll-MAY-81 
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